/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		test_P52.cpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-02-20 21:23
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-07-20 16:20
*
*   Description:
*
=============================================================================*/
#define BOOST_TEST_MODULE test_P52

#include "../../include/fem.hpp"
#include "../../include/implicit_solver.hpp"
#include "../../include/incremental_load.hpp"
#include "../../include/model/model.hpp"
#include "../../include/solid.hpp"
#include "../../include/static_equation.hpp"
#include <boost/test/included/unit_test.hpp>

LOGGING_REGISTER()

BOOST_AUTO_TEST_CASE(test_P52) {
  /*  detail, trace, process, info, warning, debug, error }; */
  geoxfem::log::log_init(geoxfem::log::severity_level::error);

  using namespace geoxfem;
  using namespace geoxfem::model;
  using namespace geoxfem::geometry;

  using MethodPolicy   = fem;
  using CouplePolicy   = solid;
  using DimensionType  = plane_strain;
  using EquationPolicy = static_equation;
  using SolvePolicy    = implicit_solver;
  using LoadPolicy     = incremental_load;
  using Model =
    geoxfem::model::Model<MethodPolicy, CouplePolicy, DimensionType,
                          EquationPolicy, SolvePolicy, LoadPolicy, mpl1>;

  std::string input_file_name =
    (boost::unit_test::framework::master_test_suite().argc < 2)
      ? "test/test_P52/P52.json"
      : boost::unit_test::framework::master_test_suite().argv[1];

  /*  detail, trace, process, info, warning, debug, error }; */
  geoxfem::log::log_init(geoxfem::log::severity_level::error);

  io::DataParser<io::json> parser;
  parser.read(input_file_name);

  /* io::DataExtractor<Model> extractor; */
  /* extractor.init(parser); */

  io::DataWriter<Model, io::vtu> writer;
  writer.init(parser.getPath(), parser.getFileName());

  /* io::DataPloter                 ploter; */

  Model model;
  /* model.append(extractor); */
  model.append(&writer);

  model.init(parser);
  model.buildProblem();
  model.solveProblem();

  MatrixXR correct_ans(12, 2), ans(12, 2);
  correct_ans << 0.0000E+00, -0.8601E-05,  //
    0.0000E+00, -0.4056E-05,               //
    0.0000E+00, 0.0000E+00,                //
    -0.5472E-07, -0.3771E-05,              //
    0.1225E-05, -0.1910E-05,               //
    0.0000E+00, 0.0000E+00,                //
    0.1792E-06, 0.5860E-06,                //
    0.6280E-06, 0.1707E-06,                //
    0.0000E+00, 0.0000E+00,                //
    0.0000E+00, 0.1129E-06,                //
    0.0000E+00, 0.1066E-06,                //
    0.0000E+00, 0.0000E+00;
  for(auto& v : shapesIn(*model.getMesh(), 0))
    ans.row(v->getIndex()) << v->toVertex().getDof("x0")->getValue(),  //
      v->toVertex().getDof("x1")->getValue();
  std::cout << "ans norm:" << correct_ans.norm() << ", " << ans.norm()
            << std::endl;
  BOOST_CHECK_CLOSE(correct_ans.norm(), ans.norm(), 1);

  auto shp0 = model.getMesh()->getShape(2, 0);
  auto shp1 = model.getMesh()->getShape(2, 1);
  auto shp2 = model.getMesh()->getShape(2, 2);
  auto shp3 = model.getMesh()->getShape(2, 3);
  auto shp4 = model.getMesh()->getShape(2, 4);
  auto shp5 = model.getMesh()->getShape(2, 5);
  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp0gp1, correct_shp0gp1,  //
    shp0gp2, correct_shp0gp2,                                  //
    shp0gp3, correct_shp0gp3,                                  //
    shp0gp4, correct_shp0gp4;

  shp0gp1 = shp0->getGaussPoint(0).getSigma();
  shp0gp2 = shp0->getGaussPoint(1).getSigma();
  shp0gp3 = shp0->getGaussPoint(2).getSigma();
  shp0gp4 = shp0->getGaussPoint(3).getSigma();
  correct_shp0gp3 << -0.4299E+00, 0.1431E+00, 0,  //
    0.1431E+00, -0.1058E+01, 0,                   //
    0, 0, 0;
  correct_shp0gp2 << -0.2511E+00, 0.8629E-01, 0,  //
    0.8629E-01, -0.6411E+00, 0,                   //
    0, 0, 0;
  correct_shp0gp4 << -0.3304E+00, 0.8352E-01, 0,  //
    0.8352E-01, -0.1016E+01, 0,                   //
    0, 0, 0;
  correct_shp0gp1 << -0.1516E+00, 0.2667E-01, 0,  //
    0.2667E-01, -0.5985E+00, 0,                   //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp0gp1.norm(), shp0gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp0gp2.norm(), shp0gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp0gp3.norm(), shp0gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp0gp4.norm(), shp0gp4.norm(), 1);

  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp1gp1, correct_shp1gp1,  //
    shp1gp2, correct_shp1gp2,                                  //
    shp1gp3, correct_shp1gp3,                                  //
    shp1gp4, correct_shp1gp4;

  shp1gp1 = shp1->getGaussPoint(0).getSigma();
  shp1gp2 = shp1->getGaussPoint(1).getSigma();
  shp1gp3 = shp1->getGaussPoint(2).getSigma();
  shp1gp4 = shp1->getGaussPoint(3).getSigma();
  correct_shp1gp3 << -0.2856E+00, 0.8499E-01, 0,  //
    0.8499E-01, -0.9141E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp2 << -0.1427E+00, 0.1394E+00, 0,  //
    0.1394E+00, -0.5806E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp4 << -0.3808E+00, 0.3735E-01, 0,  //
    0.3735E-01, -0.9549E+00, 0,                   //
    0, 0, 0;
  correct_shp1gp1 << -0.2379E+00, 0.9177E-01, 0,  //
    0.9177E-01, -0.6215E+00, 0,                   //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp1gp1.norm(), shp1gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp1gp2.norm(), shp1gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp1gp3.norm(), shp1gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp1gp4.norm(), shp1gp4.norm(), 1);

  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp2gp1, correct_shp2gp1,  //
    shp2gp2, correct_shp2gp2,                                  //
    shp2gp3, correct_shp2gp3,                                  //
    shp2gp4, correct_shp2gp4;

  shp2gp1 = shp2->getGaussPoint(0).getSigma();
  shp2gp2 = shp2->getGaussPoint(1).getSigma();
  shp2gp3 = shp2->getGaussPoint(2).getSigma();
  shp2gp4 = shp2->getGaussPoint(3).getSigma();
  correct_shp2gp3 << -0.1513E+00, 0.6412E-01, 0,  //
    0.6412E-01, -0.3680E+00, 0,                   //
    0, 0, 0;
  correct_shp2gp2 << 0.2765E-03, 0.1010E+00, 0,  //
    0.1010E+00, -0.1429E-01, 0,                  //
    0, 0, 0;
  correct_shp2gp4 << -0.2159E+00, 0.1359E-01, 0,  //
    0.1359E-01, -0.3957E+00, 0,                   //
    0, 0, 0;
  correct_shp2gp1 << -0.6432E-01, 0.5050E-01, 0,  //
    0.5050E-01, -0.4197E-01, 0,                   //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp2gp1.norm(), shp2gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp2gp2.norm(), shp2gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp2gp3.norm(), shp2gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp2gp4.norm(), shp2gp4.norm(), 1);

  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp3gp1, correct_shp3gp1,  //
    shp3gp2, correct_shp3gp2,                                  //
    shp3gp3, correct_shp3gp3,                                  //
    shp3gp4, correct_shp3gp4;

  shp3gp1 = shp3->getGaussPoint(0).getSigma();
  shp3gp2 = shp3->getGaussPoint(1).getSigma();
  shp3gp3 = shp3->getGaussPoint(2).getSigma();
  shp3gp4 = shp3->getGaussPoint(3).getSigma();
  correct_shp3gp3 << -0.2331E+00, 0.1477E+00, 0,  //
    0.1477E+00, -0.4231E+00, 0,                   //
    0, 0, 0;
  correct_shp3gp2 << -0.9447E-01, 0.1211E+00, 0,  //
    0.1211E+00, -0.9964E-01, 0,                   //
    0, 0, 0;
  correct_shp3gp4 << -0.1867E+00, 0.1015E+00, 0,  //
    0.1015E+00, -0.4032E+00, 0,                   //
    0, 0, 0;
  correct_shp3gp1 << -0.4805E-01, 0.7493E-01, 0,  //
    0.7493E-01, -0.7974E-01, 0,                   //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp3gp1.norm(), shp3gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp3gp2.norm(), shp3gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp3gp3.norm(), shp3gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp3gp4.norm(), shp3gp4.norm(), 1);

  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp4gp1, correct_shp4gp1,  //
    shp4gp2, correct_shp4gp2,                                  //
    shp4gp3, correct_shp4gp3,                                  //
    shp4gp4, correct_shp4gp4;

  shp4gp1 = shp4->getGaussPoint(0).getSigma();
  shp4gp2 = shp4->getGaussPoint(1).getSigma();
  shp4gp3 = shp4->getGaussPoint(2).getSigma();
  shp4gp4 = shp4->getGaussPoint(3).getSigma();
  correct_shp4gp3 << 0.1061E-02, -0.4210E-01, -0,  //
    -0.4210E-01, 0.7274E-01, 0,                    //
    0, 0, 0;
  correct_shp4gp2 << -0.2619E-01, -0.2217E-01, -0,  //
    -0.2217E-01, 0.9144E-02, 0,                     //
    0, 0, 0;
  correct_shp4gp4 << -0.3382E-01, -0.3302E-01, -0,  //
    -0.3302E-01, 0.5779E-01, 0,                     //
    0, 0, 0;
  correct_shp4gp1 << -0.6108E-01, -0.1308E-01, -0,  //
    -0.1308E-01, -0.5806E-02, 0,                    //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp4gp1.norm(), shp4gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp4gp2.norm(), shp4gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp4gp3.norm(), shp4gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp4gp4.norm(), shp4gp4.norm(), 1);

  ////////////////////////////////////////////////////////////////
  arithmetic::Stress<DimensionType> shp5gp1, correct_shp5gp1,  //
    shp5gp2, correct_shp5gp2,                                  //
    shp5gp3, correct_shp5gp3,                                  //
    shp5gp4, correct_shp5gp4;

  shp5gp1 = shp5->getGaussPoint(0).getSigma();
  shp5gp2 = shp5->getGaussPoint(1).getSigma();
  shp5gp3 = shp5->getGaussPoint(2).getSigma();
  shp5gp4 = shp5->getGaussPoint(3).getSigma();
  correct_shp5gp3 << -0.4855E-01, 0.3616E-01, 0,  //
    0.3616E-01, 0.1373E-01, 0,                    //
    0, 0, 0;
  correct_shp5gp2 << -0.5281E-01, 0.8267E-02, 0,  //
    0.8267E-02, 0.3775E-02, 0,                    //
    0, 0, 0;
  correct_shp5gp4 << 0.2640E-03, 0.3758E-01, 0,  //
    0.3758E-01, 0.3465E-01, 0,                   //
    0, 0, 0;
  correct_shp5gp1 << -0.4001E-02, 0.9688E-02, 0,  //
    0.9688E-02, 0.2469E-01, 0,                    //
    0, 0, 0;
  BOOST_CHECK_CLOSE(correct_shp5gp1.norm(), shp5gp1.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp5gp2.norm(), shp5gp2.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp5gp3.norm(), shp5gp3.norm(), 1);
  BOOST_CHECK_CLOSE(correct_shp5gp4.norm(), shp5gp4.norm(), 1);
}
